NASA Contractor Report 195294 




Probability Techniques for Reliability Analysis 
of Composite Materials 


Robert C. Wetherhold and Anthony M. Ucci 
State University of New York 
Buffalo, New York 

(NASA-CR-195294) PROBABILITY 
TECHNIQUES FOR RELIABILITY ANALYSIS 
OF COMPOSITE MATERIALS (State 
Univ. of New York) 45 p 


N94-37633 


Unci as 


July 1994 


^f/24 0017054 


Prepared for 

Lewis Research Center 

Under Grant NAG3-862 



National Aeronautics and 
Space Administration 




PROBABILITY TECHNIQUES FOR RELIABILITY ANALYSIS OF COMPOSITE MATERIALS 


Robert C. Wetherhold and Anthony M. Ucci 
Department of Mechanical and Aerospace Engineering 
State University of New York 
Buffalo, New York 14260-4400 


SUMMARY 

Traditional design approaches for composite materials have employed deterministic criteria for failure 
analysis. New approaches are required to predict the reliability of composite structures, since strengths and 
stresses may be random variables. This report will examine and compare methods used to evaluate the 
reliability of composite laminae. The two types of methods that will be evaluated are fast probability 
integration (FPI) methods and Monte Carlo methods. In these methods, reliability is formulated as the 
probability that an explicit function of random variables is less than a given constant. Using failure criteria 
developed for composite materials, a function of design variables can be generated which defines a “failure 
surface” in probability space. A number of methods are available to evaluate the integration over the 
probability space bounded by this surface; this integration delivers the required reliability. The methods which 
will be evaluated are: the first order, second moment FPI methods; second order, second moment FPI 
methods; the simple Monte Carlo; and an advanced Monte Carlo technique which utilizes importance 
sampling. The methods are compared for accuracy, efficiency, and for the conservativism of the reliability 
estimation. The methodology involved in determining the sensitivity of the reliability estimate to the design 
variables (strength distributions) and importance factors is also presented. 


1. INTRODUCTION 
Composite Materials 

Fibrous composite materials consist of at least two constituents: a series of purposefully oriented fibers, 
surrounded by a solid matrix. Typically, the fibers act as load-carrying members while the matrix transfers 
the load between them while fixing the fibers in the desired orientation and location within the composite. 

The resulting material is both strong and stiff. In addition to improving strength and stiffness, the mixture 
improves the survivability of the material by improving toughness and protecting the usually susceptible 
fibers from environmental damage. The most common form of composites is a stack of relatively thin layers 
of fibers and matrix (laminae), each with a specific fiber orientation, bonded together to form a laminate. 

The orientations and the stacking sequence (of laminae) can be controlled in order to generate a variety of 
mechanical properties. 

Due to the relatively low density of the matrix material, composites usually have much higher strength- 
to-weight ratios and stiffness-to-weight ratios when compared to traditional structural materials (e.g., metals). 
In addition, composites are often more fatigue resistant than these materials. Unlike traditional materials, 
which have nearly equal properties regardless of the direction of measurement, composites frequently have 
anisotropic material properties. In most unidirectionally oriented fiber-reinforced composites, the tensile 
strength and modulus are a maximum in the longitudinal or axial direction of the fibers. Many other 
properties, such as thermal and moisture expansion, thermal conductivity and impact strength are also 
directionally dependent. Most laminates use a multidirectional stacking sequence to provide a desired or 
tailored set of material properties. 

Although composites offer the possibility of tailoring material properties to a particular application, 
their design requires more care than monolithic metals or ceramics. The elastic properties of composites are 
relatively well controlled and reproducible, but the strength of composites may be more widely distributed 
due to the inherent inhomogeneity of the microstructure. It is sensible then, to speak not of "safety factors," 
but of reliability, using a knowledge of the distribution of composite strengths. We shall first review the use 
of deterministic failure theories, then indicate how a "design variable" approach can be used to calculate a 
reliability value. 
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Review of Deterministic Failure Theories 


Failure analysis in a design problem consists of the comparison of the stress (due to the applied load) 
to the resistance (or strength) of the material. With the application of multiaxial stresses, a failure theory 
based on the combined stresses must be used. In traditional materials (such as aluminum or steel) the theories 
commonly used are the maximum allowable shear stress (Tresca) theory or the distortional energy (von Mises 
criterion). These theories utilize the isotropy of the materials and are based on plastic deformation as the 
mode of failure. Composites are anisotropic and do not normally exhibit this type of failure. Therefore failure 
theories developed for isotropic materials are not applicable to composite materials. Many new failure 
theories have been developed specifically for composites; these rely on a variety of failure definitions. They 
also consider the directionality and geometry of the composite. 

For application in a thin laminated structure, the strength of the lamina is characterized by five 
strength properties in the material's coordinate system: 

SLt = Longitudinal (Axial) Tensile Strength 

Sl c = Longitudinal (Axial) Compressive Strength 

Sj t = Transverse Tensile Strength 

Sj c = Transverse Compressive Strength 

S^Ts = In-plane (Interlaminar) Shear Strength 

Both the longitudinal and transverse strength are sign-dependent, while the shear property is independent of 
the sign of the loading, given the usual orthotropic symmetry. These five properties must be determined 
experimentally. 

For failure analysis the plane stress condition of an orthotropic lamina is assumed. This is shown in 
figure 1-1. 



Figure 1-1: General Stress State of an Orthotropic Lamina [10] 
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The lamina contains unidirectional fibers at a fiber orientation angle of © with respect to the primary 
geometric axis (x). The proposed failure theories require the transformation of the applied stresses from the 
lamina or geometric coordinate system (x, y) to the principal material coordinate system (1, 2). This 
transformation is defined by equations (1.1) to (1.3). 

Cf 1 1 = c xx cos 2 0 + OyySin 2 © + 2x X ySin0cos© (1.1) 

CT 22 = Oxx s i n2 © + Oyycos 2 © -2T xy sin0cos0 ( 1 .2) 

T 12 = ( a yy “ °xx )sin©cos© + x X y (cos 2 © - sin 2 ©) (1.3) 

where a ji, 022 and x 12 are the transformed stresses, longitudinal, transverse and shear respectively, and o xx , 

<r and Xw are the applied stresses. The two failures theories used to evaluate reliability in this report are 
yy *y rr 

the Tsai-Hill theory and the more general Tsai-Wu theory. 

Tsai-Hill: This theory, developed by V.D. Azzi and S.W. Tsai, and also referred to as the Azzi-Tsai- 
Hill failure theory, is based on the anisotropic yield criterion for metals (developed by Hill). It proposes that 
failure occurs when equation (1.4) is satisfied (ref. 10): 



2 

g ll g 22 °22 

2 + 2 
c z 

^Lt *Tt 



(1-4) 


This equality is for tensile loads only. If compressive loads are applied (either in the longitudinal or 
transverse direction) then the corresponding compressive strengths are substituted. This criterion defines a 
failure envelope which is the graphical representation of the failure theory in the stress coordinate system. 
The envelope forms the boundary between the the safe region and the failure region. Using o j j and a 22 as 
the variables defining the coordinate axes the failure theory can be defined on a quadrant-by-quadrant basis. 

In the +0 1 j/ -KJ 22 quadrant, 



g ll g 22 


Lt 




In the +c 1 ] / -Ct22 quadrant. 



g H g 22 





In the -ctj j / +022 quadrant, 


(1.5) 


( 1 . 6 ) 
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(1.7) 


In the -O] 1/ -022 quadrant. 
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(18) 


2 2 2 

CT 11 _ g ll g 22 g 22 X 12 

2 2 + 2 - 1 2 
S Lc 4c Sj c S LTs 


These four equations define a failure envelope. The continuous but non-smooth nature of the failure envelope 
is due to the anisotropic strength characteristics. 

Tsai-Wu: This failure theory (developed by S.W. Tsai and E.M. Wu) proposes that failure occurs when 
the following equality (equation (1.9)) is satisfied (ref. 10): 


2 2 2 

F l a ll +F 2 0 22 + F 6 x 12 + F ll a ll +F 22 a 22 + F 66 T 12 + 2F 12 a ll a 22 = 1 


(1.9) 


where Fj and F- are called the strength coefficients and are defined as: 


F, = 


1 S Lt S Lc 


Fo = 


2 S Tt S Tc 


f 6 = 0 


F,, = 


11 ' S Lt S Lc 


F 22 S Tt S Tc 


F 66= 2 


LTs 


with Sj^ c and Sj Q > 0. The Fj 2 term represents the interaction of the longitudinal and transverse stresses. 

Since all other factors are known, this term can be determined by a number of experimental procedures. One 
possible technique is to perform an equi-biaxial test (with applied transverse and longitudinal stresses only), 
where On = O 22 = O at failure. In this case equation (1.9) reduces to: 


(F | + F 2 ) o + (F 1 1 + F 22 + 2Fj2) ^ ^ (1*10) 

Using this equation and the definitions for the strength coefficients the strength interaction term Fj 2 can be 
recovered: 


F 


12 “ 



1 - 


(t-t 


+ 






+ 


1 

S Tt S Tc 



OH) 


In order to ensure a closed failure surface, the following inequality must be satisfied: 

-\V\\*22> m ^\2*\Vu*22> m 

Due to the limited amount of experimental data available, the lower limit of this range will be used in this 
report, consistent with an analogy to the two-dimensional Von Mises criterion (ref. 14). These equations 
define the boundary between safe and failure regions. Unlike the envelope defined by the Tsai-Hill failure 
criterion, this failure envelope is continuous and smooth. This is due to the presence of linear terms in the 
failure equation. 

The envelopes defined by these two failure theories are compared graphically in figures 1-2 for a 
typical carbon fiber/epoxy matrix composite. 
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Figure 1-2: Tsai-Wu and Tsai-Hill Failure Theories[10] 


Note that the scale of the a 22 axis has been exa gg erated t0 show the curvature. The failure envelope defined 
by the Tsai-Wu failure theory is a continuous and smooth ellipse while the Tsai-Hill envelope is smooth 
within each quadrant, and only continuous at the interface between quadrants. Both envelopes intercept the 

a 1 1 axis at S^ t and — Sj^ c , and the 022 at Sj t and -S-p c . 

Application of these theories is limited to a deterministic analysis. Given a known set of stresses and 
strengths, these failure criteria can predict failure or survival. For the reliability analysis involved in a design 
problem, this is of little use. In most cases we are attempting to analyze the reliability of specimens of an 
unknown, but statistically describable, strength. The only true method to determine the strength for a given 
specimen would be to load that specimen to failure. This would give us an accurate measure of strength but 
in so doing would destroy the specimen. The reliability method desired would use deterministic stresses (we 
assume these can be accurately measured or calculated) with composite strengths following known dis- 
tributions. These composite strengths would be obtained by testing specimens until the distributions of their 
strengths could be predicted with a degree of confidence. The reliability methods discussed in this report will 
use this information and the failure theories available to predict reliability. 


Design Space/Integration Method 

To fully understand the methods used for calculating reliability, the term “design space” must be 
defined. This n-dimensional space (hyperspace), where n is the number of design variables, represents the 
entire domain of possible values of the variables. A joint density function, which predicts the relative 
probability of occurrence of a specific set of values for the design variables (coordinates in design space) is 
defined over this entire domain. The next step is to determine which areas of the design space, representing 
possible sets of design variables, will result in the failure of the system. Although we usually consider the 
strengths as design variables and the applied stresses as known values, it is straight-forward to add the 
stresses to the set of design variables. 

In general, a failure criterion defines a multi-dimensional surface which divides the hyperspace into 
two domains; a safe region and a failure region. To determine the reliability, the defined joint density 
function is integrated over the safe region while conversely the probability of failure is the integral over the 
failure region. The integral of a density function is termed a (cumulative) distribution function. For a single 
variable system this analysis requires the integration of the single-variable density function over a defined 
range. Problems with multiple variables require the more complex (multivariate) joint density function. 
Consider the miltivariate case with n=2. If the variables x and y are both normally distributed, then their 
joint density function decays continuously from the mean values (ji x ,|Ay) and appears as a “pile of sand.” 
Figure 1-3 schematically shows this for two design variables (where a x = a y = a, the standard deviation). 
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Figure 1-3: 2D Joint Density Function 


For the failure criteria defined in the previous section, the design problem consists of three variables 
for the Tsai-Hill failure criterion and five variables for the Tsai-Wu criterion. The joint density functions 
associated with these variables are similar to the one shown in figure 1-3, but generalized to three or five 
dimensions. The first step in the integration of these density functions is to define the bounding surface for 
the safe region. Once this is accomplished, the density function is then integrated over the volume (area) 
contained within the surface. The complexity of the typical bounding surface defined by these failure criteria 
makes these integrals very difficult to evaluate in closed form. Fast Probability Integration (FPI) methods 
rely on approximating the failure surface by a predetermined geometric form. These geometric forms are 
chosen specifically to match analytical forms for which evaluation of the integral for the reliability is 
practical. In choosing an appropriate approximating geometric form, a compromise between accuracy and 
computational difficulty is made. Monte Carlo methods are also widely used, although they may require 
many simulation trials using a random sampling method to evaluate the integral. 

In this report we will consider the application of both FPI and Monte Carlo methods to the calculation 
of reliability for composite materials under multiaxial loading. The composite materials considered are: a 
unidirectional fiber, polymer matrix composite (PMC) of carbon (graphite)/epoxy, and a ceramic matrix 
composite of silicon carbide (SiC)/reaction bonded silicon nitride (RBSN). The analysis will demonstrate 
the sensitivity of the calculated reliability to parameters of the strength distributions. Both normal and non- 
normal parameter distributions will be evaluated. In addition to this, a general comparison of the methods 
will be performed. 


2. RELIABILITY ANALYSIS TECHNIQUES 

The process involved in performing an integration to evaluate a probability of failure can be viewed in 
its simplest form as a design problem involving a single random variable. If the system involved has a 
deterministic resistance (strength), x, and a random applied load, Y, then the reliability of the system is the 
probability that the resistance is greater then the load, P(x > Y). Conversely the failure probability is the 
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probability that the load is greater then the strength, P (Y > x). To determine this probability the distribution 
of the load must be determined. First, define a probability density function by: 

fy(y) dy = P [y < Y < y+dy] 

Once the probability density function for the load, fy(y), is known, an integration is performed over the 
unsafe region (where the load is greater than the strength) to determine the probability of failure. This is 
shown graphically in figure 2-1: 



In the specific case of a Gaussian random variable the load Y can be “reduced” or transformed into a 
standard normal variable (SNV), y, with a mean of zero and a standard deviation of one unit: 

Y = Hy + yOy (2.1a) 

or 

y = (Y — |Xy )/a y (2.ib) 

where |iy is the mean value of Y and G y is the standard deviation of Y. This variable can then be used in the 
previous reliability equation: 


Reliability = P (x > [ly + y G y ) (2.2a) 

or 

Reliability = P ((x - p,y)/Gy > y) (2.2b) 

At this stage the variable on the right hand side of the inequality is a SNV and all the quantities on the left 
are known. These known quantities can then be reduced to a single value defined as the “reliability 
coefficient” or “reliability index,” 6. The reliability of the system can then be stated as: 
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Reliability = P (8 > y) 


(2.3) 


In the single variable case, this inequality defines a “safe region.” This region is the set of values of y for 
which the structure will not fail. The next step is to perform the integration implied by equation (2.3) over 
the appropriate region. For this case the reliability of the system is equal to the probability that a SNV is less 
than 6. This is easily determined from tables of the SNV Distribution or through calculations of any other 
known distribution function. As before the probability of failure, P y, is the complement of the reliability: 

Py = 1 - Reliability = 1 - P (8 > y) (2.4) 

Note that the use of the words “reliability index” for P usually implies that the value x exceeds the mean of 
Y (|iy), P > 0, implying a reliability greater then 0.5. The equations (2.3) and (2.4) are, however, valid in any 
circumstance. 

This process is similar when the system involved has a random resistance (strength), X, in addition to 
a random applied load, Y. The probability of failure, Py , of the system is still the probability that the load is 
greater then the resistance, P(Y > X). To determine this probability, the density and distribution of the load, 
fy(y) and Fy(y), and of the resistance, f x (x) and F x (x), must be determined. The distribution function F is 

defined by: 

F Y (y) = P [Y < y] 

y 

= j fy(z)dz 


Once these probability functions are known, an integration is performed over the unsafe region (where the 
load is greater than the strength) to determine the probability of failure. This is shown in figure 2.2: 



II Overlap/Failure m X Of Y 

y Region x 


Figure 2-2: Reliability Analysis for Two Variable System 


The required probabilities can be formulated in the following manner: 


Py = P(X < Y) = £P(X < YIY = y) P(Y = y) (2.5) 

all y 
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For continuous variables the summation becomes: 


P/= J F x (y)f y (y)dy 
o J 


An alternative expression is given by: 


P/= J [1- F y (x)]f x (x)dx 
o J 


(2.6a) 


(2.6b) 


2.1 Design Variables 

The two important characteristics that must be known about a design variable are how its values are 
distributed and if it is statistically independent of the other design variables. The variable distributions fall 
into two categories: normally and non-normally distributed. The majority of the techniques described in this 
report are for normal variables only. Non-normal variables can be used directly in the simple Monte Carlo 
method or approximated for use in the Rackwitz-Fiessler method; these special techniques will be briefly 
considered. All integration methods in this report are limited to uncorrelated variables. With known 
correlation factors between the variables, an orthogonal transformation can be performed on the variables 
such that the new variables are uncorrelated. This transformation requires the following steps for normally 
distributed variables. First the covariance matrix, [V], is determined. The elements of the covariance matrix 
are defined as: 


Vij=E[(Xi-w)(Xj-jij)] (2.7) 

E[ # ] = expected value of [•] 

where the entries represent the dependence (correlation) of one design variable on the other. Since [V] is 
positive definite, there exists an orthogonal matrix, [T] such that [T][V][T]^ is a diagonal matrix. These 
diagonal elements are the eigenvalues of the covariance matrix. We can use this transformation matrix on 
the variables to produce uncorrelated variables by the following method: 

Y = TX (2.8) 

With these new variables the covariance matrix is diagonal and therefore the variables Y are uncorrelated. 

The simulation methods can also be performed with correlated variables if the correlation is known. In our 
application (composite materials) the material strength design variables may be correlated, but the nature of 
this correlation is unknown. We shall assume that the variables of strength and applied load are uncorrelated. 


2.2 Determining The Most Probable Point of Failure 


The general case consists of a design problem involving multiple random variables, X (a vector of 
dimension n). Note that X may include both applied loads and strength (or resistance) values. To determine 
the probability of failure, a general criterion, g (X), which is a function of the design variables, is required. 
This criterion defines a surface between the “safe region,” g (X) < 0, and the “failure region,” g (2D > 0, and 
is known as the “failure surface” g (X) = 0. To facilitate calculation, the variables are converted into 
standard normal variables (SNV) by the following change of variables: 
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Xj — Xj CTj + (-t j 


i=l to n 


(2.9a) 


or 


*i = (X j - Hi )/Oj i=l to n (2.%) 

where |ij is the mean value and is the standard deviation of the random variable Xj, the ith design variable. 
The failure surface is also transformed into a function of these reduced variables, giving g(x) = 0. The 
variables, Xj, are normalized in units of standard deviation and measure the distance (in probability space) 
from the mean. The shortest distance from the origin (mean values) to the failure surface is then determined 
using a minimization technique. This distance is the reliability coefficient, 6. The closest approach point on 
the failure surface, x*, is defined to be the most probable point (MPP) of failure since the density function for 
the SNV’s decreases monotonically with distance from the origin. Using this information, the reliability can 
be calculated by employing a probability integration method which can provide a quick and accurate 
solution. 


2.3 Importance Factors 

Once the MPP has been determined, the importance of considering the variability (randomness) of an 
individual design variable can be determined. The importance is measured by the amount a variable's value 
at the MPP, Xj*, deviates from its mean, p.j , in units of standard deviation. With reduced variables, this is 
simply the absolute value of the magnitude lxj*l. To compare the relative importance of the randomness in 
the individual design variables, “importance factors”, which are simply the direction cosines, are defined: 

Ij = (Xj*/p) i = 1 to n (2.10) 


This ratio Ij of the variable’s value at the MPP to the reliability coefficient is always between one and 
negative one. A variable whose randomness does not affect the design problem (has a low importance) will 
not vary from its mean value; Xj* ~ 0. The more the importance factor deviates from this, the mean value, 
the greater the importance of the variability. These importance factors can be used to simplify a design 
problem. If a variable consistently has a factor showing low importance, then it may be possible to replace 
this random variable with its mean value in the failure criterion. This would simplify the problem by reducing 
its dimension. This is of minor value to the problems being analyzed in this report due to the small number 
(3 or 5) of variables being employed, but can be important in problems of larger dimension (n). 


2.4 Fast Probability Integration Techniques 

Fast Probability Integration (FPI) methods use a variety of techniques to evaluate the measure of 
probability space enclosed or bounded by a failure surface. This is accomplished by first determining the 
most probable point (MPP) of failure in the reduced variable probability space. The reliability coefficient, 6, 
is the distance from this point to the origin in units of standard deviation. For First Order methods, only this 
information is required. From the knowledge of 6 an approximation of the failure surface is made (a boundary 
shape is assumed) and the integration performed. For ease of calculation only simple shapes are assumed. In 
Second Order methods, additional information about the curvature of the failure surface at the MPP is 
required. These Second Order methods use both the MPP and the curvature to provide more accurate 
approximations. 


2.4.1 First Order Reliability Techniques 

Design Problems with Linear (Planar) Failure Surfaces . — In some situations the failure surface is 
planar (linearly dependent on the design variables) and an exact solution can be calculated. In these 
situations the equation for the failure surface has the following form: 
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g(X) = a jXj + a 2 X 2 ... + - y 


(2.11a) 


or 


n 

g(X)= X a i x i“Y (2-1 lb) 

i=l 

where (Xj are constant coefficients and y is a constant. The equation of the failure surface after the 
normalization of variables becomes: 


n n 

Z<Wi =Y-Z a i^i 

i=l i=l 


and the distance to the MPP, B, can be determined directly from the geometry: 


B 


n 

X a iPi- Y 
i=J 




i=l 


( 2 . 12 ) 


(2.13) 


For the applications being evaluated in this report a linear failure surface does not arise and is of minimal 
interest. In the more general case the same steps are performed, but a more robust technique is required. 

The Hasofer-Lind (H-L) Method. — This is the simplest First Order, Second Moment reliability 
method used to evaluate probability. This technique approximates the failure surface (a hypersurface) by a 
hyperplane tangent to the failure surface at the MPP. This is shown in Figure 2-3 for a system of two design 
variables in SNV space: 



Planar 

Approximation 


^ Failure Region 


X 


2 


t 

Figure 2-3: Hasofer-Lind FPI Technique 


The plane can be defined by the first derivatives of the surface (at the MPP) and the MPP in terms of 
reduced variables: 
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n 


X ( x i " x i*> (dgOO /5x i)* = 0 (2.14) 

i=l 


First, however, the MPP point must be found. The usual analytical technique is to formulate a minimization 
problem with the distance to the MPP (in units of standard deviation) as the function to be minimized: 


Minimize: 


1/2 


D(x) = 



(2.15) 


Subject to: g(x) = 0 


Note that the minimization problem 2.15 is formulated in terms of reduced variables. Once the problem is 
formulated a minimization routine will solve for the point that offers the smallest value of D()0, satisfying 
g(x) = 0. The minimization solution for D(x) gives the distance to the MPP, (6), as well as the importance 
factors (Ii). The tangent plane defined by 2.14 may now be found. However, the 6 value can be used directly 
to solve for the probability that the structure is within the safe region. By assuming the failure surface is 
planar and using reduced (SNV) variables, the value of the density function integrated over the hypervolume 
is found to be equal to the standard normal integral (distribution function) at B: 

Reliability = <J> (B) (2.16) 

while the probability of failure is the complement: 


P/ = 1-<D(B) = <D(-B) (2.17) 

where <I> is the standard normal distribution function. The evaluation of <l> either by computer or by 
comparison with tables will determine the reliability or the probability of failure. Note that equations (2.16) 
and (2.17), commensurate with figure 2-3, assume that reliability exceeds 50 percent; that is, that the origin 
is not part of the failure region. If this is not the case, a similar equation may be used. 

Hypersphere Centered at the Origin Approximation. — The biggest potential problem with the H-L 
method is in the planar approximation. In the case of a failure surface that is convex about the MPP (shown 
in fig. 2-3), the surface is termed 'friendly’ and the linear approximation is to some extent conservative. The 
alternative is a concave or ’dangerous’ curvature. This results in a liberal approximation and may 
underestimate the actual probability of failure. This situation, for a system of two variables, is shown in 
figure 2-4: 
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g(x)=0 


Planar 

Approximation 


Failure Region 


Figure 2-4: Hasofer-Lind with dangerous curvature 


One response to this problem is to approximate the failure surface by a hypersphere centered at the origin. 
This method is also categorized as first order in that it is based on the distance to the MPP only. This 
approximation is shown in figure 2-5: 



Hypersphere 

Approximation 

Failure Region 


Figure 2-5: Origin Centered Hypersphere Approximation 


The hypersphere is defined solely by the distance to the MPP, 8: 


£ x*j 2 - B 2 = 0 


(2.18) 
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This approximation is extremely conservative in that it encloses the smallest possible volume in the safe 
region, based on the calculated value of 6. By assuming the failure surface is spherical and by using SNV’s, 

the sum of the squares (X x*^) has a chi-squared distribution with n degrees of freedom. Therefore the value 
of the integration over the hypervolume can be obtained from the chi-squared distribution with n degrees of 
freedom: 


Reliability = X n 2 (B 2 ) (2.19) 

while the probability of failure is the converse: 

Pf = 1 - x n 2 (B 2 ) (2.20) 

This probability can be determined by computer or by the use of tables. 

Rackwitz-Fiessler Method. — When the design variables are not normal variables, the standard 
normalization definition {x = (X-p)/o} may still be used. There is no assurance, however, that the location 
of the closest point to the origin will be the most probable point. This difficulty arises from the fact that the 
decay in probability density may be unequal along the different design variable axes. If we can arrange for 
the probability description at the closest approach point to resemble (be locally equivalent to) a normal 
distribution, we may be able to force the closest point to be the MPP. The reliability would then be 
evaluated from the normal distribution. This is the rationale behind the Rackwitz-Fiessler method (ref. 12), 
which requires an iterative solution as outlined below. 

Let the design variables Xi have any non-normal distribution. On the initial iteration, instead of 
normalizing the variables in the usual fashion, convert the minimization problem into a function of the 
original variables. 


Minimize: 


D(X) 


X)= 


i=l 


1/2 


( 2 . 21 ) 


Subject to: g(X) = 0 

where jij and Oj are the means and standard deviations of the non-normal variables (assumed known) and 
g(X) is the failure criterion in the original variables. This problem is mathematically equivalent to the 
minimization problem mentioned previously. The only difference is that the normalization of the variables is 
done within this minimization problem while in the previous problem the variables were normalized 
externally. The next step is to match the density and distribution functions of a similar normal distribution. A 
new distribution, with mean and standard deviation (oi'), is obtained for each variable, Xj, by matching 
the original distribution and density function at the closest approach point, X* (which is equivalent to x* in 
the normalized coordinate system), with those of a normal variable; 

<D ((Xi* - Pi y<y\') = F (Xi*) (2.22) 

♦ ((Xj*- pi’)/Oi') = f (Xi*) (2.23) 


where 0(.) is the standard normal distribution function and <))(.) is the standard normal density function. To 
determine the new parameters (means and standard deviations) of the equivalent normal distributions 
equations (2.22) and (2.23) must be inverted to yield the following equalities. 


Oi’ =<t> 1 [F(Xj*)] }/f(Xj*) 

(2.24) 

Hi' = Xi* - <jj' 4> - 1 [F(Xj*)J 

(2.25) 
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The new values of the means and standard deviations found by equations (2.24) and (2.25) are inserted in the 
minimization problem 2.21, and the closest point X* is found. These two steps: solving the minimization 
problem and determining equivalent normal distributions, are repeated until a convergence criterion (either 
absolute or relative) is satisfied. The distance calculated is then used in the planar or hyperspherical 
approximation of the failure surface outlined above. From these approximations the reliability (or probability 
of failure) can be determined. 

2.4.2 Second Order Reliability Techniques 

As in the first order methods, the second order methods approximate the hypervolume of the failure 
region by a simpler volume over which the probability integration may be more easily evaluated. This is 
equivalent to determining a simpler bounding surface which approximates the true surface. To improve the 
approximation of the failure surface beyond the level employed in the first order methods, additional 
information about the failure surface is required. In the second order methods the approximation relies on the 
curvature of the surface at the most probable point of failure in addition to the distance to that point. With 
this additional information, a more correct approximation surface can be used with confidence. The only 
limit to the shape chosen is the ability to match it to a known closed form expression for probability. 

The General Quadratic Form Solution . — For implementation of any of the second order methods, the 
failure surface must be twice differentiable in the neighborhood of the checking point in the normalized 
coordinate system. Once the failure surface has been converted to the normalized coordinate system, the 
next step is to expand the failure surface, g(x) = 0, into a second order Taylor series about the checking 
point, x* (ref. 5): 


r n 


g(x) = g(x*) + 


X dgte ), (x ._ x .* ) + l (x ._ x .*)2 

3xj |x=x* ( i * ) + 21 / , a x .2 |x=x* ^ x i x i > 


i=l 


n-l n 


*1 1 


L i=l 


3 2 5(x) 


“ .4* 3x,-3bc.- \x = x 

i=i j=i + 1 1 J 




= 0 


(2.26a) 


or in matrix notation after rearrangements: 

g(x) = (X - X*) T - [G x ] • (x - X*) + 2 {g x ) T -(x - X*) + 2-g(x*) = 0 (2.26b) 

where [G x ] is a matrix of second and mixed derivatives and {g x } is a vector of first derivatives with the 
following form: 


[G x ] 


“ d 2 g 00 d 2 g00 

dx, 2 9 x 1 3x 2 

d 2 g(x) 8 2 g(2l) 
^ x 2 dx l 3 x 2 2 


~i 

3x l 5x n 


3 2 g(x) 3 2 g( 2 L) 

L^ x n 9x l 3x„ 2 J 


(2.27) 
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(2.28) 


{g x ) 


1 9g(x) 3g(x) 

1 ^ x i dx 2 



T 


By a process of linear transformations equation (2.26b) can be converted to one of two standard forms in a 
new coordinate system (defined as the z-system). These two possible transformed quadratic forms are 
defined by equations (2.29) and (2.30). 

n 

IXi (Zi-Sj^ci (2.29) 

i=l 

m n 

X *-i ( z i “ Si ) 2 + S ( c i “ z i) = c 2 (2-30) 

i= 1 i=m+l 


where Xj are eigenvalues of the [G x ] matrix, the 8 j terms are the noncentralities in the z-system, and the 
constants Cl, C 2 and cj (i=l,n) are dependent on the linear transformation. The form these equations define 
is dependent upon the signs of the eigen value-to-constant ratios, Xj/C] and Xj/C 2 » for equations (2.29) and 
(2.30) respectively. If all ratios are greater than zero, the quadratic is positive definite. In the case of 
equation (2.29) this defines an ellipsoid centered at ( 81 , 62 ,.., 8 n ) with semi-axes [(Cj/Xi) 172 , (C]/X 2 ) 1/2 ,..., 
(Ci/X n ) 172 ]. If any of the eigenvalues (or ratios) are zero, the equations represent cylindrical forms. In the 
case of eigenvalues with differing signs, the quadratic is indefinite and can have a variety of forms. 

Because the variables in the z-system are SNV's, new variables can be defined 

n 

W = X^i^i-Si ) 2 (2-31) 

i=l 

m n 

V=XM z i-*i) 2 + X(ci-Zi) (2-32) 

i=l i=m+l 

where the variable W is a combination of noncentral chi-squared distributed variables only, while V contains 
a linear combination of normally distributed variables. The distribution function of W is the probablistic 
content of a normal distribution where the spherical region is defined by equation (2.29). The distribution of 
variable V can also be defined in a similar fashion by equation (2.30), although the shape of the region is 
dependent upon the constants C 2 and cj (i=l,n). From the definition of the safe domain (in either the initial 
x-system or in the transformed z-system) it can be shown that the probability of failure is estimated by 

Pf = P(W > Cl) = 1 - F W (C]) (2.33) 

Pf = P(V > C 2 ) = 1 - F v (C 2 ) (2-34) 

where F w and F v are the probability distribution functions of W and V respectively. A variety of expressions 
have been derived for use in estimating these density functions. These approximations require complicated 
integrations and as a result have limited applicability in the study of reliability. Two simpler forms for the 
quadratic approximation have probability distributions that are relatively easy to calculate and will be 
evaluated in detail in the following sections. 

Special Forms with Predetermined Principal Axes. — In addition to their use in the exact curvature 
solution, the curvatures at the checking point or MPP x* can be used to approximate the surface with a 
predetermined shape. This results in a substantial savings in computation time but decreases the accuracy of 
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the results. The consequence of these approximations can only be evaluated when compared to the ‘exact' 
solution. 

The first step in performing this method is to determine the MPP of failure as previously described. 

The next step is to calculate the curvatures of the failure surface at this point, with an eye towards 
employing a single resultant cuvature. This is accomplished by using the direction of the principal axis, 
which is a unit vector from the origin to the MPP. The coordinate system (x-system) is rotated in such a way 
as to allow the principal axis (nth axis) to lie along the position vector from the origin to the MPP. The most 
probable point is converted from (xj*, X2*, ... , x n *X in the old coordinate system, to (0, 0, ... , 8), in the 
new coordinate system. This new system is defined as the y-system. This orthogonal rotation requires the 
transformation matrix, T, where: 


X = TY (2.35) 

This is shown in figure 2-6: 



Since we have all of the information about the failure surface in the unrotated coordinate system, the inverse 
of this transformation is of greater interest: 


Y = T- 1 X = T T X (2.36) 

The inverse and transpose are equal since T is an orthogonal matrix composed of unit vectors. To 
determine this matrix we use the MPP whose position is known in both coordinate systems. From this the nth 

row of the T^ matrix can be determined: 


T^ = Xi*/B (2.37) 

Since the (n- 1) secondary axes in the Y coordinate system are arbitrary, but must be perpendicular to each 
other and lie in the hyperplane normal to the nth axis, the rows of the matrix have similar requirements. 
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The rows must represent unit vectors which are perpendicular to each other. The first step is to solve for row 
(n-1). The first n-1 positions (in the n-1 row) are assumed to be equal to one. Then, using the perpen- 
dicularity requirement, the n^ 1 element of row (n- 1) is solved for by taking the dot product of row (n-1) with 
row (n) and setting the result to zero. This requires solving one linear equation with one unknown (element n 
of row (n-1)). The n- 1 vector is then normalized. The computed vector is both a unit vector and perpen- 
dicular to the nth row. This process is then repeated for the next lowest row, (n-2). In this case the last two 
positions (n- 1 and n) are non-arbitrary and must be solved for by taking dot products with the n- 1 and nth 
rows. This requires solving two linear equations with two unknowns. These two steps are repeated for each 
row until the entire matrix is complete. 

The second order and mixed derivatives at the MPP are calculated in the rotated coordinate system. 
This can be accomplished in two ways. The first is to use the matrix to transform the variables x into the new 
coordinate system y, then transform the failure criterion using these variables, and finally take the derivatives 
in the new system. The second method is to use the transformation matrix to determine the values of the 
derivatives in the rotated (y)coordinate system from those calculated in the unrotated (x) coordinate system. 
This requires the following calculations: 



3 2 g(x) 

fyfyly* 


9 2 g(x) t T t T 
ax p 3x lx* >'p j'q 


(2.38) 


_ dg&) _ dgQO t T 
g yi' 9ypy* 9x lx* i'P 


(2.39) 


where [G y ] is the matrix of second and mixed derivatives, {g y } is the vector of first derivatives in the rotated 
y-system and y* = T^x*. 

The coordinate system then undergoes an additional rotation with the y n = z n axis fixed, such that the 
mixed derivatives go to zero in the z coordiante system. This is accomplished by performing an orthogonal 
transformation like the coordinate rotation described in Section 2.1. This final coordinate system is defined 
as the z-system, and its curvatures are the “principal” curvatures. The first and second derivatives ({ g z } and 

[G z ]) are then calculated in the z-system using either of the two methods mentioned above. Using this matrix 
and vector, the curvature of the surface with respect to the variables in the final rotated coordinate system (z) 
can be calculated. For an n-dimensional surface there are (n-1) principal curvatures, Kj, calculated. This is 
accomplished by solving for the roots of the characteristic equation in the y-system: 

det [Gy] - K • ij = 0 (2.40) 


where gy is the first derivative with respect to y n , solution for K will give the vector of curvatures, [I] is the 

(n-1) rank identity matrix and all derivatives are performed at the MPP. In the z-system the curvatures can 
be determined directly from the ratio of derivatives since the [G z ] matrix is diagonal: 


K i = 


G z .. 


i =1 to n- 1 


NO SUM 


(2.41) 


These curvatures are identical to those calculated with equation (2.40). It is important to note that although 
use of the y- and z-systems simplify the calculation of the curvatures, these curvatures can not be identified 
with the original design variables due to the rotations performed. In both coordinate systems (y and z) a 
“dangerous” curvature (concave towards the origin) is defined by a negative second derivative while a “safe” 
curvature (concave away from the origin) is defined by a second derivative with a positive value. The first 
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derivative with respect to the primary (or nth) axis, gy , is negative by definition, since it defines the 

boundary between the “safe” region (where the failure criterion is positive) and the “failure” region (where 
the failure criterion is negative). Taking these two conditions into account the sign of the curvature is 
opposite to the sign of the second derivative used to calculate it; a “dangerous” curvature is positive and a 
“safe” curvature is negative. The radius of curvature, Rj, is the reciprocal of the curvature: 


Rj = 1/Kj i =1 to n-1 (2.42) 


Non-Central Hypersphere Approximation. — The first geometric form to be utilized as an 
approximation is the non-central hypersphere. Unlike the hypersphere used in the first order method, this 
hypersphere is not centered at the origin. This is shown in figure 2-7: 



Figure 2-7: Noncentral Hypersphere Approximation 


To determine the reliability one of the n- 1 radii of curvature values must be chosen, since a hypersphere only 
uses one radius. The choice of radius depends on the nature of the approximation desired. For this 
approximation only the positive and finite radii are considered; this eliminates any curvatures which are 
negative ("safe") or zero (locally linear). If a conservative result is desired, the smallest radius, and thus the 
smallest hypersphere, is chosen. The opposite is true if a liberal result is desired. In addition to these 
choices, the average of the curvatures can be used. Once the radius is chosen, the equation of the 
hypersphere in the z-system becomes: 

n-1 

R 2 = [z n - (6- R)] 2 + ^z^ (2.43) 

i=l 
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This hypersphere has radius R and center at the point (0, 0, ... ,6 - R), in the z-system. If the z variables are 
standard normal variables, the equation (2.43) is noncentrally chi- squared distributed with a noncentrality 
parameter 5 = [R - B]^. This is used for reliability analysis: 

Reliability = * 2 g (R 2 ) (2.44) 

and 

p /= ' - * is « 2 > < 2 - 45 > 

This distribution can be calculated from a known distribution equation (reference 11). 

Rotational Parabaloid Approximation. — This is the other quadratic geometric form used in this 
report to approximate the hypervolume contained by the failure surface. This is shown in figure 2-8: 



Figure 2-8: Rotational Parabaloid Approximation 

The same choice of an appropriate radius of curvature must be made - either the minimum, maximum or 
average. As before, this choice is made from the set of positive and finite radii. Once the radius is chosen, 
the equation of the rotational parabaloid, in the z-system, becomes: 

n-1 

4 X z f + (Zn " 6) = 0 (2 - 4 

i=l 

Since the z variables are SNV's, the first term in the equation is chi-squared distributed while the second 
term is normally distributed. Summing these two terms, the probability of failure becomes the probability of 
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a chi- squared variable convoluted with a noncentral normal variable’s being greater than zero. This results in 
the following equation: 


oo 



0 


(2.47) 


in which f x 2 is the density of a chi- squared variable with v degrees of freedom and <J> is the SNV distribution 
v 

function. In this case the density function has only n-1 degrees of freedom due to the form of the previous 
equation. Evaluating the failure probability equation (2.47) requires numerical quadrature of the specified 
distribution and density functions. 

The two presented quadratic approximations use equations based on those found in the article by 
Fiessler, Neumann, and Rackwitz (ref. 5). The specific form of the equations varies significantly from those 
referenced, however. To duplicate the results found in this report the modified equations (2.43), (2.46) & 
(2.47), as defined above, were required. The form of these equations in the article (ref. 5) are as follows: 

n-1 

R 2 = X \ + [z n " ( R + B)]2 (2 - 43 ' } 

i=l 


n-1 


2R 




(z n - B) = 0 


i=l 


(2.46') 



(2.47') 


In addition to these equations the noncentrality parameter was defined as 8 = [R + B] , and the center of the 
hypersphere was located at (0, 0, ...,0, R + B). The majority of discrepancies regarded changes in sign. 

These differences may be caused by a difference in sign convention or in the definition of the normalized 
variables. All attempts to use these primed equations as stated in the two approximations (Non-Central 
Hypersphere and Rotational Paraboloid) gave incorrect results. 


2.5 Monte Carlo Methods 

Simple Monte Carlo Technique. — This is the simplest technique available to analyze a design 
problem and calculate Ff . It consists of sampling a number of random points, where the distribution of 
sample points is controlled by the distribution of the design variables. In each iteration, a value is generated 
for each design variable. These values are then tested in the failure criterion, g(x)- The results of this test 
are consistent with our previous definition: 

g(x) < 0: Point falls within failure domain 
g(x) = 0: Point falls on the failure surface 
g(x) > 0: Point falls within safe domain 
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This is repeated, with the results of each simulation being recorded. Once a sufficient number of iterations is 
completed, the probability of failure is calculated in the following manner: 


Number of Failure Simulations _ Nf 
Total Number of Simulations N 


(2-48) 


This is a very simple approach with only one serious drawback. As the reliability increases, the probability of 
failure decreases, and the total number of simulations required for an accurate estimate of Py increases. For 

engineering design problems it is not unusual to have acceptable Pf values to be as low as lx 10“^ or 

lxl0“9. p or reasonable accuracy, a minimum of 20 to 30 sample failures are generally required, with a 
commensurate requirement on the total number of samples of say (20/P f). This is extremely costly in terms 
of computation time. 

Efficient Sampling Technique. — By using a portion of the FPI methods outlined above, a more 
efficient sampling technique has been developed to address the inherent inefficiency of the basic Monte 
Carlo technique. In an attempt to reduce the number of samples, this technique is designed to generate 
variables in a limited sampling volume that follows the original distribution (within that limited sampling 
volume). As in the FPI technique outlined above, the first step is to determine the MPP and the associated 8. 
Since this point is defined as the closest approach of the failure surface and the reduced variables are all 
SNV's, it follows that any point generated which lies at a distance (from the origin) less then 8 would fall 
within the safe domain. This n-dimensional sphere, of radius 8 and centered at the origin, termed the 6- 
sphere, can then be considered a portion of the safe domain and removed from the sampling domain. 

Removal of this domain gives a substantial decrease in the number of samples required, since most points 
generated would otherwise fall near the origin within the 8-sphere. This is shown in figure 2-9: 



Sampling Domain 


The probability content contained within the 6-sphere is calculated by using the formula for the cumulative 
chi-squared distribution with n degrees of freedom. By excluding the 6-sphere from the sampling domain a 
substantial improvement in the efficiency is obtained. 

To explain the steps used in this probability of failure estimator, the sampling domain must be split 
into the two portions defined by Ixl < 8 and IxJ > 8. The probability of failure can then be shown in the 
following equation: 

P/=Pr(g(x)<0) (2.49) 
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(2.50) 


Consider the two domains using the Theorem of Total Probability: 

P/= P{g(i) < 0 | Ixl < 6} • P{lxl < B) + P{g(x) < 0 | Ixl > 6} • P{W > B} 

Instead of using the conditions based on Ixl and 6, we can express them equivalently using Ixl 2 and B 2 for use 
with the chi- squared distribution: 

P/= P(g(x) < 0 I Ixl 2 < B 2 } ■ P{lxl 2 < B 2 } + P{g(x)<0 I Ixl 2 > B 2 } • P{lxl 2 > B 2 } (2.51) 

Replacing the various probability terms with their known distributions and values results in: 

P/= 0 • (B 2 ) + P{g(x) < 0 I Ixl 2 > B 2 } • (1- (B 2 )) (2-52) 

In the first term the first part is replaced by zero due to lack of possible failures within the B-sphere while the 
second term is the probability the point is within B-sphere. In the second term the second part is the 
probability that a point lies outside of the B-sphere. The next step is to determine the truncated probability 
density function, f( r (x), given that Ixl > B: 



Kr r -f(x) 

0 


if Ixl > B 
if Ixl < B 


(2.53) 


K tr is a constant used to scale (normalize) the distribution, defined as 

K /r = 1/(1- X„(B 2 )) ( 2 - 54 > 

such that the integral of the density function over the applicable domain equals one: 


oo 


1 f tr (x) dx = 1 
Ixl > B 


(2-55) 


The final probability term can then be replaced by this truncated density function: 


P{g(x) < 0 1 Ixl > 6} = {ffi-(x)dx (2.56) 

gOO < 0 

= Jl(x)frr(x)dx (2.57) 

all x 


where I(x) denotes the indicator variable defined as 


I 



1 

0 


if x is in the failure region 
if x is in the safe region 


(2.58) 
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The probability of failure becomes: 


P{g(x) < 0 I Ixl > 6} = E frr [ I(x) ] 


(2.59) 


where Ef^ [ I(x) ] denotes the expected value of I(x) with respect to the fj r - distribution. By combining 
equations (2.52) and (2.59) the probability of failure is: 


P/=d-Xn(B 2 ))E ffr [1(20] 


(2.60) 


Like the simple Monte Carlo method, the estimated or inferred value of the indicator variable is: 


N 


1 Nf 

Ef, r [!(*)] = ^ £*<*> = n 


(2.61) 


where N is the total number of simulations and Nf is the number of failure simulations. These simulations are 
generated according to the (truncated) f/ r - distribution. This is accomplished by converting the basic 
variables (Xi, X 2 ,...,X n ) to polar coordinates (R, ©i, ©2>— ,®n- 1) (ref. 9), where (0], ©2,.M©n-l) defines 
the direction of X while R defines the length of X in x-space. If the Xj's are standard normal variables it can 
be shown that the distance R is independent of the direction (defined by (©i, ©2>***>©n-l)) (ref, 9 pg. 18.) 
These simulations are realized by sampling a random direction unit vector a = (aj, <X2,.*.,oc n _]). This 
requires the sampling of n-1 dummy SNV’s, Y\> Y 2 ,..., Y n _i, which are then normalized to a unit vector f; in 
the following manner 

£i = jyj i=l to (n-1) (2.62) 


where IYI is the length of the vector formed by the n- 1 normal variables 


I Yl = 





1/2 


(2.63) 


This method is used to generate the (n-1) variables 0i which define the direction. The length variable, R, is 
generated following a truncated chi-squared distribution. This can be accomplished by using realizations of a 
chi-squared random variable and then discarding the values which do not fall within the reduced sampling 
domain (R > P). Due to the rapid decline of the chi- squared density function (see figure 2-10A) a majority 
of the simulations would be discarded, which is highly inefficient. An alternate method proposed by Harbitz 
(ref. 6) involves the following variable transformation: 

U = exp — R2^j (2.64) 

where a is a constant greater than two and dependent upon p. The sampling domain for U, [u\ y U2], 
corresponding to the sampling domain for R, [ri , r2l, where r2 = (ri + 3) (beyond n + 3 the chi-squared 
density function is assumed relatively insignificant), is given by: 
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[ui, U2] = 


exp [ - — ro2 
1 a 


exp 


(-H 


The density function for u is then proportional to g(u), where: 

g C0 )-[-ta(u)](“«)- 1 uta / 2)-l. 


(2.65) 


(2.66) 


This is a more uniform density function (see figure 2-10B) and therefore offers a more efficient sampling 
method. 




Figure 2-10: Sampling Density Function 

In practice the following two procedures are considered for generating a variable which follows the 
distribution of equation (2.66) and shown in figure 2-10B: 

Procedure 1: For large values of P or a small number of variables n, generate two stochastically 
independent rectangularly distributed (uniform density function) variables: 

U = Rect [ui, u2] 

G = Rect [0, g(ujxiax)] ( 2 - 67 ) 

where 


u max - ex P 



( 2 . 68 ) 


Sample random values of U and G until G < g(U). This “acceptance sampling” is required to generate 
variables which follow the desired distribution. The corresponding R-value used in the simulation is then: 


R=W-alnU 


(2.69) 


25 



This procedure is repeated until the desired number of simulations is generated based on the required number 
of sample failures. 


Procedure 2: For small values of [5 and large n (number of variables), generate R by one of the 
following two equations depending upon if the number of variables is odd or even: 

\ n ' 2 1 

For even n: R = - 2 ln| Jl^ij (2.70) 


f(n-l)/2| 

For odd n: R = - 2 In j H Vj > + X 2 (2.7 1 ) 


where Vj are independent unit rectangular variables and X is a standard normal random variable. This 
procedure is repeated until the desired number of simulations is generated. 

* 

If we now return to the reliability equation (2.60), a probability of failure estimate, , is given by: 

* 2 o Nf 

p f = 0- X n (B 2 )) • j* (2.72) 

When comparing this method to the simple Monte Carlo method, we see a substantial increase in efficiency. 
This can be shown by the following comparison of the total number of simulations required: for the simple 
Monte Carlo method N * 10/P y (assuming that 10 failed simulations gives an accurate estimate) while for 

2 o 

the Importance Sampling method N * 10 - Py/(1 - X n (B ))■ This means that for a design problem with 

Py= l.OxlCT^ and 6 = 3.5, N =100 for the Improved Sampling compared to N * 1x10^ for the simple Monte 
Carlo method. In contrast to the FPI methods proposed earlier, this method does not require any prior 
knowledge of the shape of the failure surface. 


3. RESULTS & DISCUSSION 

There are three areas of reliability analysis being evaluated in this chapter. These will be covered in 
Sections 3.1, 3.2 and 3.3. As a standard multiaxial stress state, we shall consider the loading of an "off-axis" 
specimen with orientation angle, ©. The stress state is given by a reduction of equations (1.1), (1.2) and 
(1.3): 


<*11 

r cos 2 © 

<*22 

> = \ sin 2 © 

. x 12 . 

L-sin©cos0. 


(3.1) 


The first portion of Section 3.1 will compare results obtained by using the techniques discussed in Chapter 2 
for design problems with normally distributed strengths. Accuracy of the individual methods will be 
evaluated by comparing calculated reliability estimates. These comparisons will be performed over a range 
of orientation angles and applied stresses. In addition to accuracy, the other areas of comparison will include 
computational speed, intrinsic conservatism of the methods, and degree of difficulty in implementation. This 
comparison will be designed to demonstrate the methodology involved, and as such will not be exhaustive in 
its presentation. The second portion of Section 3.1 will evaluate the Rackwitz-Fiessler method. This is the 
basic method available to analyze systems with non-normally distributed design variables. This evaluation 
will include all areas of comparison discussed in the first portion. Section 3.2 will evaluate the sensitivity of 
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the probability of failure estimates to uncertainties in the original strength distribution parameters. The 
sensitivity study will be performed in a variety of cases with changes in parameters, orientation angles and 
applied stresses. Due to the difficulty in determining a closed form solution (with the notable exceptions of 

laminae with 0° and 90° orientation angles) a numerical sensitivity study will be employed. As in the 
previous section, this analysis will be more concerned with demonstrating the methodology than in providing 
exhaustive detail. The final section (3.3) will provide and discuss the importance factors found during the 
reliability analysis. 

Confirmation of Published Solutions. — To test the computational techniques used in performing 
these methods, a series of published design problems was duplicated. For the Hasofer-Lind method, the 
problem duplicated was analyzing composite failure using Hashin’s Criterion with the applied stresses as the 
random variables [ref. 2]. These results of the published analyses gave relatively low reliability (on the order 
of 75%) and were not compared to results calculated with other methods. The results found from our analysis 
did not match the published results. Consultation with the principal author (Elishakoff) eventually led to the 
determination of the errors in these results. These errors were caused by a change of units which took place 
between the original problem [ref. 3] and the journal article form. The original problem was posed in units of 
MSI which were converted in the published problem to GPa (with a conversion factor of 6.895). This 
conversion was not performed correctly on the standard deviation values (in [ref. 2 ]: standard deviation of all 
stress distributions = 1 GPa, correct value: 6.895 GPa). Once this problem was corrected, the published 
results and our calculated results matched nearly exactly. 

Additional analysis of other design problems were also duplicated (refs. 1, 4 and 5). For the Rackwitz- 
Fiessler method the problem duplicated involved flexure of a steel beam; see Ang & Tang (ref. 1). The three 
random variables had lognormal or asymptotic extreme-value distributions. Our analysis of this design 
problem offered results similar to the published results; however, more iterations were required. This is likely 
due to Ang & Tang’s use of an analysis technique which was tailored to design problems having lognormal 
and asymptotical distributed variables, while our approach allowed for analysis of design problems containing 
variables with any nonnormal distribution. 

The second order methods (both the Noncentral Hypersphere and Rotational Parabaloid) were tested 
using analysis of a bar under a series of tensile loads (ref. 5). This design problem was analyzed by all of the 
FPI methods and the results were compared with “exact” values calculated with the Monte Carlo method. 
When this design problem was analyzed using second order methods, extremely accurate estimates were 
possible. This was attributed to two reasons. Because this was only a two dimensional design problem (two 
variables), there was only a single radius of curvature involved; thus the radius used in the analysis was the 
exact radius at the MPP. Also the radius of curvature was of the same order as the reliability coefficient. 

This resulted in a failure surface which had a relatively extreme dangerous curvature and limited the surface 
to shapes which are close to the approximations. The Importance Sampling Method was tested using a 
design problem with a failure surface which, in reduced variables, was a rotational parabaloid (ref. 6). This 
example offered both calculated probabilities of failure and information on the efficiency of the method. The 
probability of failure results analyzed with our computational techniques were close to the published results. 
The deviations between the two sets of results are likely due to the small number of simulations used (N = 
100). The portion of this example involving efficiency did not match well with the published results, 
however. It was determined that the matrix of a-values found in table 2 of reference 6 were transposed (i.e., 
the a-values found in the first row of the table actually belong in the first column). Analysis of these 
problems using the corrected a-values offered efficiency results which matched those published. 

3.1 Comparison of Techniques 

Accuracy. — The details of the comparison of lamina strength results using various techniques can be 
seen in tables 3-1 to 3-3. Each table consists of the reliability analysis of a lamina with a specific orientation 
under a series of axially applied stresses. At each stress level the reliability is analyzed by the basic Monte 
Carlo and Efficient Sampling Techniques, as well as four FPI methods: the Hasofer-Lind Method (Planar 
Approximation), the Central Hypersphere Approximation, the Non-Central Hypersphere Approximation and 
the Rotational Parabaloid Approximation. These results are compared to the reliability estimate calculated 
by the Monte Carlo method, which is believed to be the most accurate, provided sufficient trials can be 
realized. The number of simulations performed by the Monte Carlo method is chosen to give an extremely 
accurate (defined “exact”) estimate. This is done by requiring the number of simulated failures to be larger 
than 50; this requires 50/Pf simulations. 
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The material used in these simulations is a unidirectional fiber, polymer matrix composite (PMC) 
consisting of carbon (graphite) fibers in an epoxy matrix. The nominal or mean values for the five strength 
properties based on typical values are as follows: 

Longitudinal Tensile Strength: 1400 MPa 
Longitudinal Compressive Strength: 1200 MPa 
Transverse Tensile Strength: 70 MPa 
Transverse Compressive Strength: 140 MPa 
In-Plane Shear Strength: 70 MPa 

The distribution of these properties will be assumed to follow a normal distribution with a coefficient of 
variation (standard deviation/mean) of 0.10. For the second order approximations (the Non Central 
Hypersphere and the Rotational Parabaloid Approximations) a series of radii are calculated. The infinite 
radii are eliminated. The second order methods chosen are incapable of analyzing safe curvatures (Ki < 0) 
so these are also discarded. From the remaining curvatures the largest curvature (smallest radius of 
curvature) is chosen. Typically there exists only one positive and finite (defined “dangerous”) curvature for 
cases analyzed by the the Tsai-Hill and Tsai-Wu failure criteria. One zero value curvature (R = ©©) was 
usually discarded for Tsai-Hill cases, while the set of discarded curvatures for Tsai-Wu cases usually 
consisted of one zero value and two “safe” curvatures (Kj < 0). 


Table 3-1: Reliability Analysis of a 30° Gr/Ep Composite 


Failure Criterion: 

Tsai - Hill 1 

Tsai - Wu 

Applied Stress 
(MPa) 

70 

90 

120 

70 

90 

120 

MPP-Normalized 







(x*) 







SLt 

-0.0010 

-0.0016 

0.0002 

-0.0565 

-0.0705 

-0.0403 

SLc 




0.0690 

0.0909 

0.0578 

Sxt 

-0.173 

-0.310 

-0.363 

-0.455 

-0.630 

-0.390 

STc 




0.140 

0.165 

0.0906 

SLTs 

-5.52 

-4.09 

-1.70 

-5.28 

-3.69 

-1.12 

Reliability 
Index p 

5.52 

4.10 

1.74 

5.30 

3.74 

1.19 

Radius of 
Curvature R 

1190 

570 

225 

440 

235 

115 

Probability of 
Failure: 

Ctl. Hypersph.: 

108.e-8 

75.8e-5 

38.9e-2 

3.52e-5 

1.55e-2 

9.23e-l 

Planar 

Approximation 

1.67e-8 

2.02e-5 

4.13e-2 

5.83e-8 

9.04e-5 

1.17e-l 


1.78e-8 

2.29e-5 

4.96e-2 

7.57e-8 

13.7e-5 

1.63e-l 

Rotational 

Parabaloid 

1.67e-8 

2.02e-5 

4.13e-2 

5.97e-8 

9.35e-5 

1.21e-l 

“Exact” 

Value 


2.02e-5 

4. 1 3e-2 


9.73e-5 

1.33e-l 

H 1 

1.61e-8 

2.10e-5 

4.72e-2 

5.59e-8 

9.77e-5 

1.27e-l 


From table 3-1 it can be seen that the various FPI methods, with the exception of the central 
hypersphere, predict approximately the same probability of failure. This is due to the large radii of curvature 
found in all cases, which force all the models to behave like a linear (planar) approximation. It is only at the 
highest stress levels, where the radius of curvature is at its lowest value, that these predictions differ 
significantly. In all cases the central hypersphere approximation differs substantially from the “exact” value 
and thus is of interest only insofar as it offers a lower limit on reliability. In all six of the cases shown in 
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table 3-1 (three loads with two criteria) only one curvature (per case) was positive and finite or in reliability 
terms “dangerous”. This curvature was chosen for the second order evaluation. As mentioned previously, by 
choosing the largest curvature (and thus the smallest radius), the approximation becomes conservative. 
Including this curvature appears to have a slight effect on the Tsai-Wu analysis while the Tsai-Hill analysis is 
virtually unchanged. In the Tsai-Hill analysis only one curvature is ignored; using the Tsai-Wu analysis, 
three curvatures are ignored since they are either zero or negative. These extra curvatures can describe a 
surface that is significantly different than the one defined by the single “dangerous” curvature, and caution 
must be used. The results calculated by the Importance Sampling method do not depend upon curvature and 
thus are unaffected by this concern. The Non-Central Hypersphere approximation like the central 
hypersphere, tends to be conservative, but is much more accurate. Due to the flatness of the surface at the 
MPP, the probability of failure calculated by the two second order methods and the probability found from the 
planar approximation are approximately equal. 


Table 3-2: Reliability Analysis of a 45° Gr/Ep Composite 


Failure Criterion: 


Tsai - Hill 



Tsai - Wu 

J 

Applied Stress 
(MPa) 

55 

70 

85 

55 

70 

85 

MPP-Normalized 







(X*) 







SLt 

0.0067 

0.0016 

0.0004 

-0.0379 

-0.0451 

-0.0197 

SLc 




0.0727 

0.0739 

0.0446 

Sti 

-0.540 

-2.92 

-1.43 

-5.36 

-2.05 

-0.909 

Stc 




0.0835 

0.202 

0.104 

SLTs 

-5.68 

-2.93 

-1.40 

-1.07 

-2.75 

-0.953 

Reliability 
Index (3 

5.71 

4.14 

2.00 

5.47 

3.44 

1.32 

Radius of 
Curvature R 

523 

529 

360 

540 

310 

215 

Probability of 
Failure: 
Ctl. Hypersph: 

3.97e-7 

65.6e-5 

262e-3 

2.93e-5 

374e-4 

88.3e-2 

Planar 

5.75e-9 

1.72e-5 

2.28e-3 

4.73e-8 

2.94e-4 

9.31e-2 i 

approximation 







Non Central 
Hypersphere 

7.28e-9 

2.73e-2 

4.95e-2 

7.09e-8 

4.01e-2 

19.2e-2 

Rotational 

Parabaloid 

5.75e-9 

2.26e-4 

4.54e-2 

2.29e-8 

3.01e-4 

9.46e-2 

“Exact” 

Value 


5.13e-5 

3.03e-2 


4.27e-4 

10.6e-2 

Importance 

sampling 

1.29e-8 

6.21e-5 

3.52e-2 

5.01e-8 

3.14e-4 

9.81e-2 


In table 3-2 the results calculated by the various FPI methods (except the Central Hypersphere) are 
extremely similar. Only the failure probability values calculated by the Planar Approximation are 
systematically below the “exact” values. This is expected because this approximation is liberal in dangerous 
curvature cases. All other reliability methods (second order FPI and Importance sampling) appear to have 
equal accuracy. This is significant in the first two simulations performed using the Tsai-Hill Failure Criterion 
(Applied Axial Stress of 55 and 70 MPa respectively) in that the radius chosen for the approximation was the 
smaller of the two dangerous curvatures calculated. Therefore even a conservative second order 
approximation may result in a relatively accurate reliability analysis. It is important to note the change in 
the nature of the curvature as a function of the applied stress level using the Tsai-Hill failure Criterion. This 
shows that shape of the failure surface can change based solely on the stress applied. Because assumptions 
about this shape may be used when determining which reliability method should be chosen, this property 
must be considered. 
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Table 3-3: Reliability Analysis of a 60° Gr/Ep Composite 


Failure Criterion: 


Tsai - Hill 

1 


Tsai - Wu 

i 

Applied Stress 
(MPa) 

55 

60 

65 

55 

60 

65 

Normalized 







MPP (x*> 







SLt 

0.00052 

0.00057 

0.00058 

-0.0141 

-0.0134 

-0.0115 

$Lc 




0.0273 

0.0263 

0.0231 

Sit 

-3.70 

-3.03 

-2.33 

-3.49 

-2.77 

-2.02 

STc 




0.0588 

0.0657 

0.0649 

SLTs 

-0.342 

-0.385 

-0.396 

-0.53 

-0.561 

-0.523 

Reliability 
Index P 

3.72 

3.06 

2.37 

3.53 

2.82 

2.09 

Radius of 
Curvature R 

1100 

909 

766 

907 

755 

640 

Probability of 
Failure: 







Ctl. Hypersph: 

3.13e-3 

2.51e-2 

1.33e-l 

2.82e-2 

1.58e-l 

4.99e-l 

Planar 







Approximation 

9.96e-5 

1.12e-3 

8.99e-3 

2.05e-4 

2.38e-3 

1.84e-2 

Non Central 
Hypersphere 

11.4e-5 

1.32e-3 

10.7e-3 

3.20e-4 

3.87e-3 

2.95e-2 

Rotational 

Parabaloid 

11.4e-5 

1.32e-3 

10.7e-3 

3.20e-4 

3.85e-3 

2.93e-2 

“Exact” 

Value 

10.9e-5 

1.25e-3 

10.9e-3 

2.31e-4 

2.75e-3 

2.28e-2 

Importance 

Sampling 

12.06-5 

1.43e-3 

1 1.8e-3 

1.16e-4 

1.42e-3 

1.24e-2 


It can be seen from table 3-3 that the nature of the curvature is strongly dependent upon the orientation 
angle for the Tsai-Hill failure criterion. It has been found by a variety of simulations that any Tsai-Hill 
reliability analysis performed with an angle of orientation greater then 50° will result in all curvatures being 
safe. This is apparent in table 3-3. In this type of case the probability of failure determined by the planar 
approximation would be conservative. The only method capable of providing a potentially more accurate 
result is the the Importance Sampling method. However, in the three simulations of this type performed the 
Importance Sampling method did not offer an improvement in the accuracy of the calculated results over FPI. 
Nearly all of the demonstrated reliability methods performed equally well when the Tsai-Wu failure criterion 
was used. 

The one type of situation in which the Importance Sampling Technique does not perform well is that of 
low reliability (6 approaches zero). As B decreases the volume contained within the p -sphere also decreases 
and the efficiency of this method approaches that of the basic Monte Carlo method. This may be of little 
importance in application of this method to actual design problems in that low reliability cases will be rarely 
modelled. 

Two additional orientations where these techniques can be compared is 0° and 90°. Unlike the 
previous cases, these orientations can result in design problems with closed form solutions. This is due to 
simplifications that these orientations cause in the failure criteria used. At these orientations the Tsai-Hill 
criterion reduces to: 


For 0°: 



a xx °yy 


Lt 



(3.2a) 
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Even though these equations are not as simple as those derived from the Tsai-Hill failure criterion, when they 
are used the reliability estimates calculated by the Hasofer-Lind method are still quite accurate (within 5% 
of the “exact” value). 

Conservatism. — All FPI methods being evaluated in this report may suffer from a predetermined 
conservatism or liberalism depending upon the nature of the design problem being analyzed. In each case the 
decision of which method to use must be made. This decision will determine if the calculated probability of 
failure will be liberal (potentially too low an estimate) or conservative (usually higher than necessary). In all 
cases the Central Hypersphere approximation offers a conservative estimate but at the same time is usually 
highly inaccurate (see tables 3-1 to 3-3). The planar approximation on the other hand can be either liberal or 
conservative depending on the curvature of the surface at the MPP. When a failure surface has all safe 
curvatures then this approximation will be conservative; while if even one curvature is dangerous this 
approximation may be liberal. The second order approximations are even more complicated due to the 
variety of possible curvatures used to approximate the surface shape. Using a large curvature (small radius) 
is conservative, while a small curvature (large radius) is liberal. These second order methods appear to offer 
the closest approximation to the actual surface however. 

Computational Speed. — The speed of the methods is difficult to describe due to their case 
dependency. Computational speed is also dependent upon the initial guess used to determine the MPP. 
Because all methods, with the exception of the basic Monte Carlo method, require the determination of (3, 
this may be ignored. The two important items which affect the speed of the FPI methods are the order of the 
approximation and the difficulty in evaluating the mathematical expression used to approximate the 
probability of failure. Once (3 has been determined, the first order methods require the evaluation of the 
normal distribution function for the planar approximation and the chi-square distribution function for the 
central hypersphere. Both of these functions are relatively simple to evaluate computationally. 

Additionally the second order methods require the determination of a rotation matrix, the evaluation of 
first and second derivatives at the MPP, and the determination of the eigenvalues of the characteristic 
equation to determine the curvatures. These steps add substantially to the computational time. The two 
mathematical functions we must evaluate for these methods are more complicated than those used for the 
first order. The Non-Central Hypersphere approximation requires the evaluation of the non-central chi-squared 
distribution function. This is done by integrating the chi-squared density function using a numerical 
quadrature where each density function value is represented by the summation of an “infinite” number of 
terms. Accurate results for this calculation are only possible for failure surfaces (g-functions) with small radii 
at the checking point; for large radii, generation of accurate results are not possible using available 
computational techniques. For the Rotational Parabaloid the calculation of the distribution function requires 
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the integration by numerical quadrature of the convolution of a normally distributed term and a chi-squared 
term. These two methods are far more computationally time consuming then either of the first order methods. 

The Monte Carlo methods available are by far the slowest computationally, due to their repetitive 
nature. This can be shown in the following example: 

Example 3-1. — As mentioned in the Introduction, the choice of the strengths as design variables is 
based on a knowledge of composite manufacturing techniques, while the choice of the stresses as 
deterministic values is relatively arbitrary. It is not difficult to imagine a design problem in which the 
applied stresses are also random variables. In this situation the problem would have six random variables if 
analyzed with the Tsai-Hill failure criterion and eight if analyzed with Tsai-Wu. Analysis of this type of 
problem with the FPI methods given would not require substantially more computational time. The only 
change would be the difficulty of the search for the MPP. The speed of the Monte Carlo methods on the 
other hand would be greatly affected. Because the number of random variables that are required during each 
simulation is directly related to the amount of computational time required to perform each simulation, this 
change would increase the amount of time required to perform the basic Monte Carlo proportionately. The 
Importance Sampling drastically cuts down the number of simulations required but at the same time requires 
the determination of p as in the FPI methods. It is still far quicker than the simple Monte Carlo in nearly all 
cases. 

Implementation Difficulty. — The simplest method to program is the basic Monte Carlo method. All 
that is required is random number generation to obtain realizations of the strengths (and stresses, if required), 
followed by an acceptability test using the failure criterion. The first order methods and the Importance 
sampling are only slightly more difficult. These require using a minimization routine with the failure criterion 
as the equality constraint. Once this is performed and the reliability coefficient, 6, is known, the next step is 
either the analysis of a distribution function for the first order methods, or a series of random numbers 
generated in the reduced domain for Importance Sampling. These two methods are still relatively simple to 
program if a packaged minimization routine is available. By far the most difficult methods to implement are 
the second order methods. This is due in part to the number of steps involved in each method and the 
complexity of the individual steps. Once the reliability coefficient is known the coordinate system must be 
rotated to align one of the new axes with the vector to the MPP. This is followed by calculation of the first 
and second-order derivatives which are used to determine the curvatures. An eigenvalue problem is solved 
for the principal curvatures. These principal curvatures are examined and a single positive curvature is 
chosen to determine the radius of curvature based upon the conservativism desired. The parameters [J and R 
are used to evaluate the distribution functions which govern the particular approximation. Another factor that 
affects the implementation difficulty is the number of program changes required to perform each case. As the 
methods get more complicated more information must be known about the problem definition. The only case- 
dependent requirement of the first order FPI and Monte Carlo methods is the failure criterion equation. The 
Importance Sampling method requires this criterion and a parameter for the truncated distribution used to 
generate random numbers. The second order FPI methods require the equations of the first and second 
derivatives of the failure criterion at the MPP in addition to the failure criterion equation for each case. 

These additional steps must be examined in the context of the computational speed and accuracy desired. 

Non-Norma! Variables. — A brief comparison of methods can also be made for design problems 
involving independent variables with non-normal distributions. The accuracy of the probability of failure 
estimates determined by the individual methods will be evaluated over a range of orientation angles and 
applied stresses. Computational speed, intrinsic conservativism, and the difficulty of implementing the 
methods will also be compared briefly. In this comparison two FPI methods will be examined; the Rackwitz- 
Fiessler method with a planar approximation and the Rackwitz-Fiessler method with the Central hypersphere 
approximation. These results will be compared to the reliability estimates calculated by the Hasofer-Lind 
method using normal distributions that are “equivalent” (possess the same means and standard deviations) to 
the original non-normal distributions and to the basic Monte Carlo method, which is believed to be the most 
accurate. The number of simulations performed by the Monte Carlo method will be chosen to give a very 
accurate (defined “exact”) estimate. As in the normal variable section this is done by requiring the number 
of failures calculated (in each simulation) to be larger than 50. 

The material used in this analysis is a ceramic matrix composite (CMC). The specific material 
chosen consists of SiC/RBSN. This CMC is evaluated by examining its three strength properties with 
nominal (average) values: 

Longitudinal Ultimate Strength: 680 MPa 

Transverse Strength: 27 MPa 

In-Plane Shear Strength: 53 MPa 
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These properties are all for tensile loading of the material (the shear strength is independent of the loading 
sign). From experiments, the distribution of these properties has been shown to follow a Weibull distribution 
with the nominal parameters: 


Shape Location 

Parameter Parameter 

Longitudinal Ultimate Strength: 5.2 741 MPa 

Transverse Strength: 10.9 28 MPa 

In-Plane Shear Strength: 7.5 56 MPa 


Due to the limited amount of experimental data for compressive strength, this composite will be evaluated 
by the Tsai-Hill failure criterion only. The results of the accuracy analysis can be seen in the following table. 


Table 3-4: Tsai-Hill Reliability Analysis of CMC Composite 



30° I 

45° 


I 60" 

1 


40 MPa 

50 MPa 

20 MPa 

30 MPa 

12 MPa 

20 MPa 

MPP: s Lt 

579 MPa 

625 MPa 

653 MPa 

653 MPa 


653 MPa 

653 MPa 

STc 

17.8 MPa 

25.7 MPa 

10.2 MPa 

15.7 MPa 


9.04 MPa 

15.2 MPa 

SLTs 

21.0 MPa 

24.8 MPa 

52.6 MPa 

51.9 MPa 


53.0 MPa 

52.8 MPa 

Reliability 
Coefficient: p 

4.09 

2.87 

4.16 

2.93 

4.45 

3.03 

Probability of 
Failure: 

Ctl. Hypersph.: 

8.02e-4 

4.08e-2 

6.06e-4 

3.54e-2 

1.88e-4 

2.73e-2 

Planar 

Approximation 

2.55e-5 

2.02e-3 

1.58e-5 

1.70e-3 

I 

1.24e-3 

“Exact” 

Value 

4.67e-4 

2.54e-3 

4.64e-4 

9.29e-3 

1.92e-4 

9.28e-2 

“Equivalent” 

3 

3.44 

2.83 

4.01 

2.68 

4.28 

2.79 

“Equivalent” 

Hasofer-Lind 

4.85e-5 

2.31e-3 

3.10e-5 

3.72e-3 

9.25e-6 

2.64e-3 


From table 3-4 it can be seen that neither of the two FPI methods available show systematic accuracy. 
The only pattern in the inaccuracy of the methods is that the planar approximation always gives probability of 
failure (Pf) results which are smaller than the “exact” value (under-estimated) while Pf calculated with the 
Central Hypersphere is always higher (over-estimated). In terms of application to actual design problems 
these calculated Pf estimates would only be valuable as limits: the central hypersphere approximation being 
the lower limit and the planar approximation the upper. The values in the final two rows of the table are 
calculated by the Hasofer-Lind method using normal distributions “equivalent” (having the same means and 
standard deviations) to the original Weibull distributions. In all cases this relatively unsophisticated 
approach, with a planar approximation, calculates probability of failure values that are closer to the “exact” 
value than those calculated with Rackwitz-Fiessler method. However, these estimates are still too liberal for 
general application. It appears from the cases examined that as the load increases (probability of failure 
increases) the planar approximation becomes more accurate. To test this hypothesis three additional cases 
were tested: an orientation angle of 30° with an applied axial load of 60 MPa, an orientation angle of 45° 
with an applied axial load of 40 MPa and an orientation angle of 60° with an applied axial load of 60 MPa. 

In each case the probability of failure estimated using the Rackwitz-Fiessler method with a planar 
approximation was about 65% of the “exact” value. This is slightly more accurate than the previous cases. 
Why these results appear to follow this pattern is unknown, however. 
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These two approximations used in the Rackwitz-Fiessler method require nearly identical amounts of 
computational time. Both require the same iterative approach to calculate the reliability coefficient. The 
only difference is the evaluation of the two distribution functions. The chi-squared distribution function 
required for the Central Hypersphere is slightly more difficult to evaluate than the normal distribution function 
used in the planar approximation. The difference is minimal, however. These approximations have the same 
intrinsic conservativism when applied to the Rackwitz-Fiessler method as the dangerous curvature case does 
with the other FPI methods. The planar approximation is liberal while the central hypersphere is conserv- 
ative. Computational techniques used to evaluate these two approximations are also equally difficult to 
implement. 


3.2 Sensitivity 

Mean Strengths. — Any variation in the mean value of the strength distribution can have a strong 
effect on the probability of failure of the system. The rate that this variation changes Pf gives the sensitivity 
of the probability of failure to the mean value. To analyze this sensitivity a series of reliability estimates are 
made by varying the mean value while all other properties (orientation angle, applied load and other 
distribution parameters) are held constant. The calculated failure probabilities are then compared to 
determine the relationship between Pf and the mean strength values. From observation it can be determined 
that the log of Pf is roughly linearly dependent upon the mean value of strength. This is to be expected in 
that all of the distribution functions used in the FPI methods are to some extent exponential in form. To show 
how this type of numerical sensitivity study could be performed the following example is provided: 

Example 3-2. — In this example we determine the sensitivity of the probability of failure to the mean 
value of the shear strength. Preliminary results have shown that varying the mean shear strength results in the 
greatest change in the probability of failure (as compared to the axial or transverse strengths). The lamina 
being analyzed is a graphite fiber/epoxy matrix lamina with a 45° orientation angle. The coefficient of 
variation is 0.10. The other parameters for the strength distributions are taken to be the nominal values listed 
in Section 3.1. The applied stress is 90 MPa. 

This case will use the Tsai-Hill failure criterion. From the previous section, it has been shown that the 
Hasofer-Lind method gives accurate results for design problems of this type and thus will be used for this 
analysis. 

A linear plot does not present a simple relationship. If the data are viewed on a semi-log plot (log 
scale for the probability of failure) then an exponential relationship can be seen; see figure 3.1. 




Mean Shear Strength (MPa) 


Figure 3-1: Sensitivity to Shear Strength Mean Value (Tsai-Hill failure criterion) 


To analyze this case a linear relationship is assumed between the log of Pf and the mean shear strength 
value. The equation of the line determined by a least squares method is 

Log (P/) = 13.99 - 0.2678 ps LTs (3.5a) 

or 

P/= io(13.99-o.2678 hs LTs } (3.5b) 


where M^SlTs * s mean va ^ ue the shear strength distribution in MPa. With this equation a relatively 

accurate prediction of the probability of failure can be made at any mean strength value as long as it remains 
close to the values used to generate equation (3.5b). 

The sensitivity Pf to changes in the mean strength value can be evaluated using the first derivative of 

equation (3.5b): 

-^i- = 0.616X 10< 13 " - 0 2678 “SLTs> (3.6) 

d ^ s LTs 

Standard Deviation of Strengths. — The probability of failure is also sensitive to changes in the 
standard deviations of the strength distributions. To analyze this sensitivity, a series of Pf calculations are 
made while changing the value of the standard deviation with all other properties (orientation angle, applied 
load and other distribution parameters) kept constant. The failure probabilities are calculated using a 
reliability method and then compared. From this comparison, a numerical relationship between the 
probabilities of failure and the standard deviation is determined. A semi-log plot is again expected to yield 
the most insight into the results. The following example is provided to show how this type of numerical 
sensitivity study would be performed. 
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Example 3-3. — In this example, we determine the sensitivity of the probability of failure to the 
standard deviation value of the shear strength distribution. The lamina being analyzed is of graphite 
fiber/epoxy matrix with a 30° orientation angle. The coefficient of variation for all unchanging distributions is 
0.10. The means of the strength distributions are the nominal values listed in Section 3.1. The applied load is 
90 MPa. 

This example uses the Tsai-Hill failure criterion. From Section 3.1 it has been shown that Importance 
Sampling is the only method which gives accurate results for this type of design problem. The relationship 
between the standard deviation and probability of failure can be seen in figure 3-2. 




Standard Deviation of Shear Strength (MPa) 


Figure 3.2: Sensitivity to Shear Strength Standard Deviation. 

Although these values do not lie quite as smoothly on a line as in figure 3-1, a linear relationship 
between the log of the probability and the standard deviations of transverse strength can still be found by a 
least squares approach. This relationship is: 


Log (P/) = -12.45 + 1.104 os LTs 


(3.9a) 


or: 


P / =10 < - 12 ' 45+ 1 104O SLTs> 


<3.9b) 


where os LTs is the standard deviation of the shear strength distribution in MPa. This gives relatively 

accurate predictions of the probability of failure in the neighborhood of the nominal values. The sensitivity of 
the probability of failure is the first derivative with respect to the standard deviation. For the shear strength 
the sensitivity is 


-^- = 2 . 539 X 10(- 12 - 45 + 1 104 ° S LT s ) 

dcSLTs 
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Applied Stress. — The probability of failure is also sensitive to changes in the deterministic input of 
the design problem. In this report the two variables that have been considered deterministic are the applied 
stress and the orientation angle. We expect the variability in the orientation angle to be quite small, since 
the lay-up process is usually accurate. On the other hand, the applied stress may have uncertainty associated 
with it. To analyze the sensitivity of the probability of failure to the applied load, a series of Pf estimates are 
made by changing the value of the applied stress with all other properties (orientation angle and strength 
distributions) held constant. The calculated failure probabilities are then compared. A numerical relationship 
between the probabilities of failure and the applied stress is determined from this comparison. A semi-log 
plot is chosen for clarity of representation. The following example is provided to show how this type of study 
would be performed. 

Example 3-4. — In this example, we will determine the sensitivity of the probability of failure to the 
axially applied stress, g xx . The lamina being analyzed is a graphite fiber/epoxy matrix with a 30° 
orientation angle. The coefficient of variation for all distributions is 0.10. The mean values of the strength 
distributions are the nominal values listed in Section 3.1. The nominal value of applied stress is 90 MPa. 

This example uses the Tsai-Hill failure criterion. The second order method using the Rotational 
Parabaloid approximation gives accurate results for this type of design problem and will be used for this 
analysis. The relationship between the applied stress and probability of failure can be seen in figure 3-3. 
Typical of brittle materials, there is a great sensitivity to the applied stress. 
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Figure 3-3:Sensitivity to Applied Load. 

This relationship is analyzed by the same approach as in the previous examples with the following results 

Log (Pfi = 0.145 Gxx - 17.68 (3.12) 

Pf = 10< 0145 «xx “ 17.68) (3.13) 

where o X x IS the applied axial load in MPa. 
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3.3 Importance Factors 


Importance factors are a measure of the impact a specific variable's randomness has on the overall 
design problem. They are also the direction cosines of the MPP position vector in reduced variable space. 
The direction cosine associated with a variable is the ratio of the reduced variable value at the MPP to the 
distance to the MPP, which distance is the reliability coefficient p. The importance factors are defined as 
follows 


iLt = Axial Tensile Strength Importance Factor 
Il c = Axial Compressive Strength Importance Factor 

Ij t = Transverse Tensile Strength Importance Factor 

lj c = Transverse Compressive Strength Importance Factor 
Ij_js = In-plane Shear Strength Importance Factor 

The results of the importance factor analysis are shown in tables 3-5 to 3-7. Each table consists of the results 
calculated from reliability analysis of a specific lamina with a given orientation angle, using both the Tsai- 
Hill and Tsai-Wu failure criteria. The strength parameters for the graphite/epoxy composite analyzed are the 
same as those used in section 3.1. The stress values applied coincide with those used to generate the 
probability of failure values found in tables 3-1 to 3-3. The importance factors are found for the three strength 
distributions used for the Tsai-Hill failure criterion and the five used for the Tsai-Wu failure criterion. 


Table 3-5: Importance Factors for a 30° Gr/Ep Composite 


Failure Criterion: 


Tsai - Hill 



Tsai - Wu 

1 

Applied Stress, 
MPa 

70 

90 

120 

70 

90 

120 

ILt 

dLc) 

1.8xl0 -4 

3.9x 10 -4 

l.lxKT 4 

l.lxHT 2 

1.3xl0- 2 

1.9x10-2 

2.4x10-2 

3.4x10-2 

4.9x10-2 

Ijt 

dTc) 

3.1xl0“ 2 

7.6xl0- 2 

21xl0- 2 

8.6x10-2 

2 . 6 x 10-2 

17x10-2 

4.4x 10 -2 

33x10-2 

7.6x10-2 

iLTs 

*1.0 

*1.0 

0.977 

0.996 

0.987 

0.941 


From Table 3.5 it can be seen that the randomness of the shear strength dominates all of the other 
strength variabilities for this lamina. This domination is nearly complete in the two cases analyzed with the 
Tsai-Hill failure criterion, with applied loads of 70 and 90 MPa. In contrast the variability of the axial 
strengths (both tensile and compressive) have nearly no effect on the results. This lack of importance can 
also be seen in the transverse tensile strength when using either failure criterion. The variability of the 
transverse compressive strength on the other hand has a moderate effect (in Tsai-Wu analyzed problems 
only) on the design problem. The meaning of these importance factors can be seen in following examples: 
Example 3-5a. — In the last case (Tsai-Wu failure criterion with applied axial load of 120 MPa) there 
are a range of importance factors. We can compare the effect of these factors by replacing a single variable 
with its mean value in the failure criterion and reducing the number of variables in the design problem by 
one. In this example, the axial tensile strength Su is replaced by its mean value (1400 MPa). The 
reliability coefficient calculated with the variable is 1.19 while the coefficient with the deterministic value is 
still 1.19. This shows that the randomness of a variable with a importance factor of 0.034 is insignificant to 
the design problem. 
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Example 3-5b. — We continue to use the last case but in this example replace the transverse tensile 
strength variable, which has an importance factor of 0.33, by its mean value (70 MPa). The reliability 
coefficient calculated with the variable is 1.19 but the coefficient with the deterministic value is now 1.25. 
This is only a 5 % change in the reliability coefficient but results in a -11% change in the probability of 
failure estimated with the Hasofer-Lind method. This shows the significance of the randomness of a variable 
with a factor slightly larger than say 0. 1 . 

Example 3-5c. — In this example, using the previous results, we replace the shear strength variable, 
which has an importance factor of 0.941, by its mean value (70 MPa). The reliability coefficient calculated 
with the deterministic value is 54.1 This is 45 fold change in the reliability coefficient and results in a 
probability of failure which cannot be easily estimated with the Hasofer-Lind method due to computer 
limitations. This shows the overwhelming dominance of the randomness of a variable whose importance 
factor is above say 0.9. 

The lack of importance of the axial strength’s variability is expected due to the dependence of the 
failure criteria upon the reciprocal of the strengths squared. The relatively large nominal values of the axial 
strength (when compared with other strengths) reduces their impact on the failure criteria. The lack of 
importance of the transverse strength randomness may be due to the small orientation angle and loading 
scheme used. 


Table 3-6: Importance Factors for a 45° Gr/Ep Composite 


Failure Criterion: 


Tsai - Hill 



Tsai - Wu 

.._l 

Applied Stress, 
MPa 

55 

70 

85 

55 

70 

85 

iLt 

(iLc) 

1.2xl0~ 3 

3.9x1 0 -4 

2.0x1 0" 4 

6.9xl0- 3 

13x10-3 

1.3xl0 -2 

2.1x10-2 

1.5x10-2 

3.4X10- 2 

iTt 

OTc) 

95xl0- 3 

0.705 

0.715 

0.980 

15x10-3 

0.596 

5.9x10-2 

0.689 

7.9xl(r 2 

iLTs 

0.995 

0.708 

0.700 

0.196 

0.799 

0.722 


Shear strength variability continues to dominate the cases found in table 3-6. Additionally the 
randomness of the transverse tensile strength does have a strong effect on the cases analyzed with the Tsai- 
Wu failure criterion. This overall increase in importance (when compared with the results from table 3.5) is 
due to the larger orientation angle. The variabilities of the remaining three strengths (axial tensile, axial 
compressive and transverse compressive) have very slight importance. The importance of randomness of both 
the axial and transverse strengths increases with applied stress while the shear strength variability importance 
decreases. 
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Table 3-7: Importance Factors for a 60° Gr/Ep Composite 


Failure Criterion: 


Tsai - Hill 

H 

Tsai - Wu 

1 

Applied Stress, 

55 

60 

65 

55 

60 

65 

MPa 







iLt 

=0 

=0 

=0 

0.004 

0.004 

0.005 

dLc) 




0.008 

0.009 

0.011 

iTt 

0.996 

0.992 

0.986 

0.988 

0.980 

0.968 

dTc) 




0.017 

0.023 

0.031 

iLTs 

0.092 

0.126 

0.167 

0.150 

0.199 

0.250 


The importance factors of the variables shown in table 3-7 vary somewhat depending upon the failure 
criterion employed for the analysis. The variability of the transverse tensile strength dominates the cases 
analyzed with both the Tsai-Hill failure criterion and the Tsai-Wu failure criterion. The variability of the 
shear strength also has a modest contribution. The increased importance of the randomness of the transverse 
strengths is due to the increase in orientation angle. The randomness of the axial strengths has nearly no 
effect on these cases. 
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4. SUMMARY AND CONCLUSION 


In this report we examined a variety of methods used to evaluate the reliability of composite laminae. 
The initial portion of this report discussed the necessary steps involved in posing the design problem in the 
required form. The specific form used the composite strengths as design variables and the Tsai-Hill and Tsai- 
Wu failure criteria for failure determination. The first method type discussed was the fast probability 
integration (FPI) methods. These consisted of two groups: first order, second moment reliability methods and 
second order, second moment reliability methods. The second type discussed was the Monte Carlo methods. 
The two methods discussed in this category were the basic Monte Carlo and the Importance Sampling 
methods. These methods were examined in detail and their procedures for implementation were outlined. 

The next portion of this report used calculated probability of failure estimates to evaluate these 
reliability methods. A general comparison of the accuracy of the probability of failure estimates, the 
efficiency of the methods at determining these estimates, and the conservativism of the estimation was 
performed. This was designed to compare how these methods performed when analyzing the specific type of 
design problem used in this report and to describe the general methodology involved in making a comparison. 
The sensitivity of the probability of failure estimate to variations in the parameters of the distribution of the 
design variables was discussed. This involved specific sensitivity studies and a description of the 
methodology required. The final section of this report evaluated importance factors. The importance factors 
of the variables chosen for the example design problems were found, and the impact of eliminating low 
importance variables was evaluated. 

The conclusions derived in this study come in three parts: those made about the application of the 
methods discussed, those about the sensitivity analysis results and procedure, and those involving importance 
factors. The sections below will discuss these topics in terms of their general application and their 
application to the specific problems outlined within this report. 

The Hasofer-Lind method is by far the most efficient method discussed. It is not only the simplest 
method to implement, but it also requires the least amount of computational time to perform. The limitation 
of this method would appear to be the relatively coarse approximation of the failure surface. From previous 
examples (see ref. 5) this appeared to make this method highly inaccurate. From these results, the 
probability of failure estimates given by this method are as accurate as any of the methods chosen (see tables 
3-1 to 3-3). This is most likely due to the typically large radius of curvature (when compared to the 
reliability coefficient) found in these results. In other applications this may not be the case. Other possible 
benefits of this method are its insensitivity to the number of variables in the design problem and its ability to 
handle complex failure criteria. Due to these properties and the simplicity of the implementation procedure, 
this method is suggested for the initial analysis of any design problem. Steps should be taken, however, to 
verify the accuracy of the calculated results. This verification procedure will likely require the 
implementation of one of the higher order FPI (or sampling) methods. 

The Central Hypersphere method is as efficient as the Hasofer-Lind method, but offers none of its other 
benefits. In all cases analyzed, probabilities of failure calculated with this method differed substantially from 
the “exact” value. This is likely due to the flatness of the failure surface in our design problem. In addition 
to the inaccuracy of this method it also suffers from a strong sensitivity to the number of variables used in the 
design problem. The positive aspects of this method are that it is both simple and computationally quick and 
will calculate results which are always conservative. The amount this method over-estimates the probability 
of failure is the real concern. In general application, this method may result in a structure which is over- 
designed. This is both inefficient and, typically, expensive. As a result this method is suggested for 
determining conservative reliability estimates only while another reliability method should be used as the 
primary method for analyzing reliability. 

The second order methods (both the Rotational Parabaloid and the Non-Central hypersphere) offer 
some benefits which are not found in first order methods. In dangerous curvature cases these methods can 
predict probability of failure values which are closer to the actual value than those predicted by either of the 
first order methods (see tables 3-1 and 3-2). This benefit may be eliminated, however, if the approximated 
failure surface has a variety of curvatures which differ in both sign and magnitude (see table 3-1, Tsai-Wu 
criterion). These methods as presented here were limited to analyzing design problems which have 
dangerous curvatures. In addition to these limitations these methods require additional knowledge of the 
failure surface and may be sensitive to the radius of curvature chosen for the surface approximation. These 
methods are also more difficult to implement and require more computational time to perform. These 
properties limit the application of these methods to specific types of design problems. However, from our 
data, it is shown that use of the non-central hypersphere approximation with the smallest radius of curvature 
offers a reliability estimate which is conservative but more accurate then results calculated using the central 
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hypersphere method. Using these methods for design problems with radii of curvature which are greater than 
one order of magnitude larger than the reliability coefficient will not offer substantial improvements over a 
first order method, because the failure surface is so “flat” that linear surface is a good approximation. In 
general, it should first be determined that the design problem does not have R »B before the additional 
steps involved in this method are performed. These two approximations (Rotational Parabaloid and the Non- 
Central hypersphere) appear to offer very similar results. From these results it does not appear necessary to 
analyze a design problem using both approximations. 

The Importance Sampling method is far less computationally efficient than any of the FPI methods 
discussed. It does, however, offer a possible improvement in the probability of failure estimates for design 
problems with safe curvatures (see table 3-1, Tsai-Wu criterion). This is not possible with FPI methods, 
because all the methods utilize “dangerous” curvatures to make approximations. Another benefit of this 
method is that it does not suffer from any intrinsic conservativism as some of the FPI methods do. The 
probability of failure estimates calculated with this methods have a similar accuracy for all types of design 
problems. The only negative feature of this method is that it depends upon repetitive simulations. This means 
that as the probability of failure decreases the amount of required computation time increases. However, if the 
number of design problems being analyzed is small, the amount of computational time required for each 
analysis may not be prohibitive. This method is suggested if it has been determined that the FPI methods are 
giving inaccurate results or if the first or second derivatives of the failure surface are difficult to obtain. 

The final method examined in this thesis is the Rackwitz-Fiessler method. This is the only method 
capable of analyzing design problems with non-normally distributed variables. This method is slightly more 
difficult to implement than the other first order methods. It does not appear to give accurate probability of 
failure estimates when compared with the "exact" value. In fact an equivalent design problem using normally 
distributed variables (which have the same means and standard deviations as the non-normal distributions) 
appears to be more accurate. The inaccuracy of the Rackwitz-Fiessler method may be due to the type of 
design problem being used in the analysis, however. The method itself has not been found to have any 
intrinsic problems which would result in erroneous results. Due to the relatively limited examination of this 
method only the most conservative approaches can be suggested. The Central Hypersphere approximation is 
recommended for calculating an upper limit to the probability of failure. 

Sensitivity. — All of the presented methods can be examined to determine the sensitivity of the 
calculated probability of failure to changes in the parameters of a design variable's distribution. This 
procedure could be used to determine the importance of obtaining a good estimate of the variable's 
distribution. The only problem with the proposed procedure is its case dependence. To determine the 
sensitivity of a design problem to a specific set of design parameters requires a complete sensitivity study. It 
is suggested that if only one parameter is to be modified, then a sensitivity study using this method is 
appropriate; but if multiple parameters are to be modified then this procedure may be too time-consuming to 
perform. 

Importance Factors. — The importance factors section of this report offers an approach to simplify 
design problems. In general the randomness of any variable which has an importance factor approximately 
equal to zero will have very little effect on the calculated results. Using this information, the problem can be 
simplified by replacing this variable with its mean value in the failure criterion. In all of the cases analyzed 
the variability of the shear strength was found to be very important while the variabilities of axial strengths 
were found to be relatively unimportant (see tables 3-5 to 3-7). The importance of the other strengths' 
variabilities was dependent upon the applied load, orientation angle and failure criterion used. A variable 
which was found to be unimportant could be replaced with a deterministic value, which simplifies the problem 
(see example 3-5a). This procedure may not offer significant computational time savings when using an FPI 
method, due to the small number of design variables used. Unless a problem could be reduced to a single 
random variable, which would offer a closed form solution, the savings would be minimal. For simulation 
(Monte Carlo) methods, however, proportional savings in computer time could be achieved for the random 
variables eliminated.. 
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